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I.  INTRODUCTION 

In  a  previous  memorandum  report  [1]  numerical  solutions  of  the  Klein-Gordon  equation 
were  developed  with  a  view  toward  solution  of  relativistic  photoionization,  or  other  relativis¬ 
tic  quantum  electronics  problems.  While  the  Klein-Gordon  equation  captures  much  of  the 
relevant  physics,  especially  for  moderately  heavy  ions  (Z  -C  137),  it  does  neglect  the  spin 
polarization  of  the  electron.  This  memo  parallels  [1],  but  replaces  the  Klein-Gordon  equation 
with  the  Dirac  equation,  which  is  the  fully  relativistic  wave  equation  for  spin  one-half. 

The  notational  and  representational  conventions  used  herein  follow  those  of  Ref.  [2] . 


II.  CYLINDRICAL  BOUND  STATES 


A.  Starting  Equations 


As  in  Ref.  [1],  cylindrical  atoms  are  of  interest  as  a  reduced  model  that  is  sometimes  the 
only  one  that  is  computationally  tractable.  The  starting  point  is  the  symmetrical  form  of 
the  Dirac  equation 

[7m  (Pm  -  QAn)  ~m\ip  =  0  (1) 

where  are  the  contravariant  Dirac  matrices,  /3/t  is  the  covariant  operator  of  four- 
momentum,  is  the  covariant  four-potential,  q  is  the  charge  of  the  particle,  and  m  is 
the  mass.  In  two  dimensions  this  expands  to 

(i7°<90  +  +  i72<92  —  q^°A0  +  q^lAl  +  qrfA2  —  rn)ip  —  0  (2) 


Here  it  must  be  remembered  superscripts  are  contravariant  vector  indices.  We  seek  station¬ 
ary  states  in  the  form 


iP  = 


(3) 


where  ip  =  tan“1(a:2/a:1),  and  £  is  a  diagonal  4x4  matrix  with  integer  entries.  Specializing 
further  to  the  case  of  a  uniform  magnetic  held  A  =  ^pB^e^  gives  the  time  independent 
equation 


7  °uj 


7  °qA° 


0up 


?-Vf 


rn  +  —T+R+  + 


WiL  )  etivijj  =  0 


(4) 


where 


R± 


dp  ± 


qByp 

2  T 


i 

p 


(5) 
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r±  =  i71±72  (6) 

Here,  p2  =  x2  +  y2.  At  this  point  it  is  convenient  to  choose  a  representation.  We  will  use 
the  standard  representation  in  which 


where  i  G  {1,2,3}  and  rt  are  the  Pauli  matrices.  Taking  advantage  of  the  fact  that  diagonal 
matrices  commute,  it  is  easy  to  show  that  all  the  terms  dependent  on  <p  cancel  if  and  only 
if  l\  —  (.2  —  @3  —  —  1-  If  these  conditions  are  satisfied,  then 

(l°u  -  7  V°  -  m  +  h+R+  +  i>  =  0  (8) 

This  is  the  eigenvalue  equation  describing  the  radial  part  of  the  wavefunction.  It  decouples 
into  two  independent  systems,  one  for  -0°  and  (the  spin-up  solution)  and  one  for  -01  and 
-02  (the  spin-down  solution). 

B.  Spin  and  Angular  Momentum 

Due  to  the  symmetry  of  the  cylindrical  atom,  the  spin-up  and  spin-down  solutions  are 
easily  transformed  into  one  another.  Multiply  the  fundamental  two-dimensional  time  de¬ 
pendent  equation  (2)  on  the  left  by  71y3.  Use  {7M,7I/}+  =  2 to  move  71y3  to  the  right 
of  the  operator  in  parenthesis.  Make  a  change  of  variables  from  x1  to  =  —x1.  Gather  the 
factors  on  the  right  as  a  new  wavefunction 

^\xo,e,x2)=iV^°,-e,x2)  (9) 

If  the  components  of  A ^  all  have  definite  parity  with  respect  to  x1 ,  ip'  satisfies  the  original 
equation,  with  possible  changes  in  the  signs  of  A A  For  A0  and  A2,  odd  parity  results  in  a 
sign  change.  For  A1,  even  parity  results  in  a  sign  change. 

Now  consider  a  cylindrical  Dirac  atom  in  a  constant  magnetic  field.  The  field  config¬ 
uration  has  A0  even,  A1  even,  and  A2  odd.  Thus  the  scalar  potential  of  the  transformed 
problem  is  unchanged,  but  the  magnetic  field  is  reversed.  The  sign  of  (f  is  also  reversed,  due 
to  the  coordinate  reversal,  but  this  is  more  conveniently  viewed  as  a  transformation  of  i. 
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Specifically,  if  there  is  a  stationary  state  ^(x0 ,  p,  tp]  —B0,£),  then  there  is  also  a  stationary 
state 


(10) 


(  V,1(p;-^o)  \ 

-B0) 

4>3(p]  -B0) 

\  -B0)  J 

where  £'0  =  —£\,  £[  =  —£o,  £'2  =  —^3,  and  £':i  =  —£2.  Note  that  a  spin-up  solution  is 
transformed  into  a  spin-down  solution  and  vice-versa. 

It  is  instructive  to  consider  the  angular  momentum  of  the  system.  The  operator  of  total 
angular  momentum  is 

jz  =  \pz  -  idv  (11) 


where  the  first  term  is  due  to  the  spin  (see  Ref.  [2]).  The  eigenvalues  are  jz  =  £0  +  \  in  the 
case  of  a  spin-up  state,  or  jz  =  £\  —  \  for  a  spin-down  state.  If  the  transformation  (10)  is 
applied,  the  sign  of  jz  is  reversed.  The  total  angular  momentum  is  a  good  quantum  number. 
The  orbital  angular  momentum  is  a  bad  quantum  number  because  the  stationary  states 
have  £0  £3  and  £\  ^  ^2)  so  that  they  are  not  eigenfunctions  of  d^.  Any  spin-up  bound 

state  is  fully  specified  by  a  radial  quantum  number,  nr  (see  below),  and  the  total  angular 
momentum  jz .  In  the  following  we  treat  only  the  spin-up  state,  since  the  spin-down  state 
can  always  be  obtained  by  the  transformation  (10).  A  specific  spin-up  state  is  denoted  by 
the  ket  | nr,jz) 


C.  Analytical  Solutions  for  Coulomb  Potentials 


In  the  case  of  a  Coulomb  potential,  A0  =  Q / p,  where  Q  is  the  effective  nuclear  charge  [11], 
it  is  possible  to  obtain  an  analytical  solution  in  the  limit  of  a  weak  magnetic  held.  Define  Q 1 
and  Q2  such  that  "0o  =  Vm  +  ojpse~xp{Qi  +  Q2)  and  =  i\Jm  —  u>pse~Xp{Qi  —  Q2),  where 
A  =  y/m2  —  ca2  and  s  =  — 1/2  +  yj j2  —  q2Q2.  Then  one  obtains  the  following  uncoupled 
equations  for  Q\  and  Q2- 


m2(jj2p3 

jzmujcp - - - h 


pQi  +  [  1  +  2cr  —  2A p  H — - 

mujcp 


miucp 


•  _  7  _  mLOcp 

J  z  2 


Q'i  + 


*  _  7  _  mujcpz 

k  J z  2 


—  2A  (cr  + 


Qi  =  0 
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mujcp 


pQ 2  +  (  1  +  2cr  —  2A p  H — 

m 2+?p3  mujcp(a  —  Zw)  /  mw/ 

j2mwcp - - - +  - - - — 3-  -  2A  |  1  +  a+Z<JJ  + 


ry  nrUx)cpA 

2 


mujcp 


Q'-2  + 
<?2  =  0 


q  I  7  _  i  I  7 

Jz  i  2  \  J  z  \  zjm  2 

where  Zm  =  qQm/ A,  =  qQui/X,  a  =  s  +  1/2,  and  u;c  =  qB0/m.  In  the  low  held  limit 

\muc\p 2  <  4|j,|  and  |mu;c|p2  <C  2\jz  ±  Zm|  assuming  p  is  bounded.  Supposing  further  that 
\a±Zu\/\jzT  Zm |  <  jA  results  in 


+  (  _  2A  )  Qi  + 

p 


nrcncj2 - (a  +  ZJ 

P 


Qi  ~  o 


(12) 


Q2  +  (  — — - 2A  )  Q2  + 


2A  . 

mucjz - (1  +  a  +  4 

P 


Q2  ~  0 


The  solutions  that  are  finite  for  p  =  0  are 

( — 1  +  2ZW)A  +  (1  +  2cr  )/3 


Qi  =  Cxe^-^F 


2(5 


,  1  +  2cr,  2(5 p 


(13) 


(14) 


Q2  =  C2e^pF 


(1  +  2ZUJ)\  +  (1  +  2  <y)(5 

2^ 


,l  +  2(x,  2(5  p 


(15) 


Ci 

c2 


a 


(16) 


jz  5~  Z.m 

where  F(a,  b,  z )  is  the  confluent  hypergeometric  function  [12]  and  (5  =  \J\2  —  mucjz.  The 
energy  levels  for  the  bound  states  are  found  by  demanding  that  the  wavefunction  should 
vanish  as  p  — >  00.  Then  the  first  argument  of  F  must  be  a  negative  integer  or  zero.  Thus, 
the  energy  eigenvalues  for  a  weak  held  B0,  total  angular  momentum  jz,  and  radial  quantum 


number  nr  are  determined  from 


(+1  +  2Z^)A  +  (1  +  2(7  )(5 

2/3 


=  —nr 


(17) 


where  nr  G  {0, 1,2, .. .}.  An  important  point  to  note  is  that  for  B0  ^  0,  it  is  not  always 
possible  to  satisfy  both  equations  simultaneously  (see  discussion  of  Zeeman  splitting  below). 
The  validity  of  the  assumptions  used  to  obtain  (17)  cannot  be  decoupled  from  the  particular 
state  considered.  Sufficiently  far  from  the  axis  (p  — >  00),  the  approximations  involving  the 
smallness  of  thujc \ p2  will  always  break  down.  It  is  expected,  therefore,  that  states  whose 
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FIG.  1:  First  few  energy  levels  of  cylindrical  (solid)  and  spherical  (dashed)  hydrogen- like  xenon 
according  to  Dirac  theory  with  Bq  =  0. 


wavefunction  decays  rapidly  with  p  are  approximated  best.  These  are  generally  the  states 
with  the  lowest  energy. 

In  the  case  where  B0  =  0  the  energy  eigenvalues  can  be  expressed  as 

r  i  -1/2 


u  —  m 


1  + 


q2Q2 


(nr  +  y/j\  ~q2Q2') 


(18) 


Because  of  the  square  root,  the  cylindrical  Dirac  ion  (with  perfect  Coulomb  potential)  can 
only  be  in  an  s-state  provided  q2Q 2  <  1/4.  In  terms  of  the  atomic  number,  Z  =  — Q/q , 
assuming  the  orbiting  charge  is  an  electron,  this  condition  is  expressed  as  Z  <  68.  This 
differs  from  the  spin  zero  (Klcin-Gordon)  case  where  no  s-state  is  possible  at  all  for  a 
cylindrical  Coulomb  ion. 

Fig.  1  compares  the  first  few  energy  levels  in  a  cylindrical  hydrogen-like  xenon  ion  ( Z  = 
54)  with  those  of  the  corresponding  spherical  ion.  The  qualitative  structure  of  the  energy 
levels  is  similar,  but  in  cylindrical  geometry  each  energy  level  is  lowered  relative  to  its 
spherical  counterpart.  The  fact  that  the  energy  is  lowered  is  fortuitous,  because,  as  will  be 
shown  below,  a  soft-core  potential  can  be  used  to  raise  the  energy  of  the  cylindrical  state  so 
that  it  becomes  commensurate  with  the  corresponding  spherical  state. 
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D.  Numerical  Solutions  for  Soft  Core  Potentials 


In  order  to  compute  the  stationary  states  of  a  Dirac  ion  in  a  magnetic  field  without  any 
approximation,  one  may  return  to  the  fundamental  Eq.  (8)  and  solve  it  numerically.  The 
primary  difficulty  is  that  the  Coulomb  potential  has  a  singularity  which  is  not  amenable  to 
discretization.  One  solution  that  is  often  employed  is  to  use  a  soft-core  potential 

A°=  -=2= 

\J  p1  +  Sp2 

where  5p  is  a  constant  called  the  soft  core  radius.  The  soft  core  potential  converges  to  a 
Coulomb  potential  as  5p  — >  0.  Depending  on  the  state  considered,  a  very  small  value  of  Sp 
might  have  to  used,  in  conjunction  with  a  similarly  small  numerical  cell  size,  in  order  to 
satisfactorily  approximate  the  Coulomb  solution. 

Let  each  component  of  i^(p)  be  discretized  on  a  sequence  of  mesh  points  {p*  =  (2 i  — 
l)4r|f  G  N}.  Evaluating  the  derivatives  in  Eq.  (8)  via  centered  finite  differences  gives  a 
41V  x  41V  block  tridiagonal  matrix  equation,  where  N  is  the  number  of  mesh  points  used. 
This  breaks  into  two  independent  2N  x  2N  systems,  one  for  spin  up,  the  other  for  spin 
down.  The  spin  up  part  is 


f  o+  (2( 

\-iT+  D.  )  )  \<PZ  ) 

where  T±  are  tridiagonal,  and  D±  are  diagonal.  The  non-zero  elements  of  T±  and  D±  are 


T± 


,i,i—  1 


l 

2Ap 


(21a) 


r-p  1/2  T  jz  ^qBoPi 

~  -  ±  - X - 

Pi  2 


(21b) 


T±,i,i+ 1 


1 

2A ~p 


(21c) 


D±,i,i  =  qA ?  ±  m  (21d) 

where  A®  =  A°(pi).  The  spin  down  part  is  obtained  from  the  same  system  with  the  signs  of 
jz  and  B0  reversed,  and  the  column  vector  replaced  by  ('ll;1,  'ip2).  The  spin  down  solutions 
obtained  from  the  transformation  (10)  are  identical,  as  exemplified  in  Fig.  2.  Standard  sparse 
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FIG.  2:  Spin  down  wavefunction  computed  directly  (solid  lines),  and  using  the  transformation  (10) 
(squares).  The  parameters  of  the  numerical  atom  are  Z  =  54,  A p  =  0.01,  5p  =  0.1,  jz  =  1/2,  and 
nr  =  0.  The  magnetic  field  is  hujc/mc 2  =  3.4. 

matrix  packages  can  solve  the  eigensystem  for  various  subsets  of  the  eigenvalue-eigenvector 
pairs.  A  similar  system  can  be  constructed  for  spherical  atoms. 

An  important  question  is  how  small  the  soft-core  radius  must  be  in  order  to  achieve 
a  given  error  tolerance  with  respect  to  the  ideal  Coulomb  solution.  The  ideal  Coulomb 
solution  is  known  exactly  in  the  case  where  Bq  =  0.  Fig.  3  shows  the  convergence  of  the 
numerical  solution  to  the  Coulomb  solution  for  the  lowest  three  bound  states  of  cylindrical 
hydrogen-like  xenon  ( Z  =  54).  The  cell  size  and  the  soft-core  radius  are  set  equal  in  all 
cases.  Not  surprisingly,  the  |0,  state  is  the  hardest  to  resolve,  requiring  5p  =  10-4  to 
reduce  the  error  to  less  than  1%.  For  comparison,  the  characteristic  size  of  the  |0,  state 
for  Z  =  54  is  1/A  «  1.2. 


E.  Zeeman  Effect 

In  the  case  of  the  Klcin-Gordon  equation  for  spin  zero  it  is  possible  to  obtain  an  analytical 
description  of  the  Zeeman  effect,  to  first  order  in  the  magnetic  held.  In  the  case  of  spin  1/2, 
the  situation  is  more  complicated.  As  noted  above,  the  approximations  used  to  obtain 
Eq.  (17)  are  subtly  dependent  on  the  particular  state  considered.  Moreover,  in  a  magnetic 
held,  there  may  be  no  choice  of  radial  quantum  number  that  simultaneously  forces  both 
hypergeometric  functions  in  the  solution  to  zero  as  p  — >  oo.  This  situation  is  alluded  to 


FIG.  3:  Convergence  of  the  three  lowest  cylindrical  bound  state  energies  of  hydrogen-like  xenon 
with  diminishing  soft  core  radius,  5p.  The  energy  in  a  perfect  Coulomb  potential  is  denoted  ujq. 

in  Refs.  [3,  4].  On  the  other  hand,  we  are  always  able  to  find  well  behaved  numerical 
solutions  in  a  soft-core  potential.  One  explanation  for  this  is  that  the  requirement  that  both 
hypergeometric  functions  vanish  as  p  — >  oo  is  a  sufficient,  but  not  necessary  condition,  on 
the  normalizability  of  the  wavefunction. 

Based  on  the  information  in  Fig.  3,  one  can  check  the  accuracy  of  Eq.  (IT)  against  the 
numerical  soft-core  solution,  given  a  sufficiently  small  soft-core  radius.  Fig.  4  shows  such  a 
comparison  for  the  |0,  ±|)  state  with  5p  =  10-5  and  for  the  1 1,  ±4)  and  1 0,  ±|)  states  with 
5p  =  10~3.  The  approximations  used  in  the  analytical  formula  work  surprisingly  well  except 
for  the  1 1,  |)  state. 


III.  TIME  DEPENDENT  SIMULATIONS 

A.  Numerical  Integration  Technique 

In  the  standard  representation  of  the  symmetrical  form  of  the  Dirac  equation,  the  “elec¬ 
tron”  components  and  ip1)  depend  only  on  spatial  derivatives  of  the  “positron”  compo¬ 
nents  {ip2  and  ip3),  and  vice-versa.  This  suggests  an  efficient  time  integration  scheme  where 
the  electron  and  positron  components  are  leap-frogged  over  each  other.  To  see  how  this 
works,  it  is  convenient  to  put  the  Dirac  equation  in  the  Hamiltonian  form  id Qip  =  Hip  (this 
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huc/mc2 


hu>c/mc2 


Pnoc/mc 2 


FIG.  4:  Weak  field  (dashed)  and  numerical  (solid)  solution  for  energy  levels  in  cylindrical  hydrogen¬ 
like  xenon  for  the  states  (a)  (b)  |l,±^),  and  (c)  |0,±|).  The  upper  branches  are  for 

positive  jz  and  the  lower  branches  are  for  negative  jz.  For  reference,  huc/mc 2  =  1  corresponds  to 
a  magnetic  field  of  44  TG. 


in  no  way  scrambles  the  components  of  the  wavefunction) .  The  Hamiltonian  is  conveniently 
expressed  in  2  x  2  block  form  as 


H  = 


(22) 
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where 

D±  =  qA°  ±  m 


(23) 


-id i  -  d2 


-idi  +  d2  id3 


-A1  +  iA2 


-A1  -  iA2  A3 


In  order  to  express  the  leapfrog  scheme  in  a  compact  and  transparent  way,  the  Hamiltonian 
is  expressed  as  H  =  H <n)  +  H{ 12)  +  H ^  +  H^22\  where 


H <n)  = 


H{  21)  = 


D:  0 


HW  = 


Hi  22)  = 


Dehne  (j)  =  and  y  —  (t/2,?/3)-  Suppose  </>  is  known  at  time  level  n  —  1/2,  while  y  is 

known  at  time  level  n.  The  wavefunction  is  advanced  in  two  steps  using 

W  +  At/2)  =  e-iH^At/2e-iH^)Ate-iH^)At/2^t)  (27) 

^(t  +  At)  =  e-iH^At/2e-iHWAte-iHWAt/2^t  +  A</2)  (28) 


Here,  ij}(t  +  T)  must  be  interpreted  with  care.  If  T / At  is  a  half  integer,  the  argument  gives 
the  time  at  which  (j)  is  known,  with  \  being  known  a  half-step  earlier.  If  T / At  is  an  integer, 
the  argument  gives  the  time  at  which  y  is  known,  with  (j)  being  known  a  half-step  earlier. 

The  integration  scheme  of  Eqs.  (27)  and  (28)  amounts  to  a  carefully  chosen  factorization 
of  the  operator  exponential,  e~lHAt.  Leapfrogging  0  and  y  corresponds  to  factorizing  the 
upper  and  lower  rows.  Each  step  in  the  leapfrog  scheme  is  further  factorized.  In  particular, 
the  off-diagonal  elements  are  split,  and  arranged  so  that  the  two  halves  bracket  the  diagonal 
element.  This  can  be  shown  to  eliminate  the  lowest  order  factorization  error  that  would 
occur  if  the  off-diagonal  elements  were  kept  together. 

The  formal  expressions  (27)  and  (28)  have  to  be  further  approximated  to  be  useful  compu¬ 
tationally.  In  this  work  the  factors  corresponding  to  the  off-diagonal  part  of  the  Hamiltonian 
are  replaced  by  the  lowest  order  Taylor  expansion.  This  gives  the  explicit  formulas 


<P(t  +  At/2)  =  exp(-iD+At)  <j>(t  -  At/2)  -  iS'y(t)—  -iSy(t)  — 


x(t  +  At)  =  exp(— iD_ At)  y(t)  —  iS(f>(t  +  At/2) —  —  iS(f>(t  +  At/2)  — 


(29) 

(30) 
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hid /me2 


FIG.  5:  Comparison  of  energy  levels  in  a  cylindrical  soft  core  potential  with  Z  =  54  and  Bo  =  0.2, 
as  determined  by  time  dependent  (solid)  and  time  independent  (dashed)  calculations.  The  broad 
features  in  the  time  dependent  result  are  only  visible  on  a  log  scale. 

This  requires  only  explicit  finite  difference  evaluations  and  two  special  function  evaluations. 
Note  that  the  diagonal  part  of  the  Hamiltonian  is  treated  such  that  unitarity  is  preserved 
to  machine  precision.  The  entire  scheme  turns  out  to  preserve  unitarity  to  several  digits  of 
precision,  as  is  demonstrated  below. 

B.  Validation  Against  Zeeman  Effect 

A  general  method  for  determining  the  energy  levels  of  a  numerical  atom  using  a  time 
dependent  code  was  given  in  Ref.  [1],  The  method  is  used  here  to  test  the  algorithm  for 
integrating  the  time  dependent  Dirac  equation  given  above.  The  energy  levels  falling  out  of 
the  time  dependent  calculation  have  to  be  generated  for  a  soft-core  potential,  and  so  are  most 
conveniently  compared  with  the  numerical  time  independent  calculation.  This  comparison 
validates  the  time  dependent  calclulation  provided  the  time  independent  calculation  itself  is 
validated.  This  latter  validation  has  already  been  provided  in  the  form  of  comparisons  with 
analytical  energy  levels  in  a  Coulomb  potential  and  constant  magnetic  field. 

In  carrying  out  the  time  dependent  calculation,  a  random  wavefunction  is  initialized  in 
a  soft  core  potential  with  Z  —  54  and  dp  =  0.1.  The  number  of  grid  cells  is  2000  x  2000, 
the  number  of  time  steps  is  220,  the  cell  size  is  0.05  x  0.05,  and  the  time  step  is  0.02.  The 
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most  important  parameter  for  this  purpose  is  the  number  of  steps,  which  needs  to  be  large 
in  order  to  resolve  narrowly  separated  spectral  lines.  As  can  be  seen  in  Fig.  5,  the  peaks 
of  the  time  dependent  spectrum  neatly  line  up  with  the  dashed  lines,  which  denote  the 
discrete  energies  predicted  by  the  time  independent  calculation.  The  vertical  axis  spans 
several  orders  of  magnitude  so  that  the  effects  of  a  finite  temporal  window  can  be  seen.  The 
relative  strength  of  each  spectral  line  is  a  function  of  the  point  in  space,  ro,  where  the  time 
dependence  of  the  wavefunction  is  analyzed.  In  fact,  if  r0  happens  to  be  near  a  null  of  a 
particular  eigenfunction,  the  associated  energy  eigenvalue  may  not  appear  in  the  numerical 
spectrum.  In  Fig.  5,  the  diagnostic  point  is  r0  =  e!  +  e2,  where  e*  are  Cartesian  basis  vectors. 

C.  Relativistic  Ionization  Example 

One  important  application  of  numerical  solution  of  the  time  dependent  Dirac  equation  is 
in  modeling  relativistic  photoionization  processes.  In  Ref.  [1]  an  example  is  given  using  the 
spin  zero  Klcin-Gordon  equation.  In  the  Klein-Gordon  case,  the  parameters  of  the  ioniza¬ 
tion  problem  are  the  radiation  frequency,  cuo,  the  peak  vector  potential,  A0,  and  the  atomic 
number  of  the  hydrogen-like  ion,  Z .  In  the  case  of  the  Dirac  equation,  an  additional  param¬ 
eter  must  be  considerd:  the  orientation  of  the  electron  spin  with  respect  to  the  radiation 
polarization  and  wavevector.  In  the  two-dimensional  geometry  considered  in  this  report,  it 
is  only  possible  to  consider  the  case  where  the  radiation  polarization  is  orthogonal  to  the 
spin.  The  form  of  the  vector  potential  is 

A(r,  t )  =  Aq  [cos  (k  •  r  —  uj0t)  —  1]  Q(u>ot  —  k  •  r)e^  (31) 

where  0  is  the  Heaviside  step  function,  and  k  •  =  0.  In  three  dimensions,  nothing 

constrains  the  choice  of  relative  to  the  spin.  In  two  dimensions,  cannot  have  a 
component  in  the  ignorable  direction,  for  otherwise  a  momentum  would  be  induced  in  that 
direction,  leading  to  a  contradiction.  Hence,  if  the  spin  is  chosen  to  be  in  the  ignorable 
direction,  the  radiation  and  electron  polarizations  must  be  orthogonal.  In  this  report,  k  ||  e1; 
e_R  =  e2,  and  the  spin  is  aligned  with  e3.  The  electrostatic  potential  is  a  soft  core  potential, 
with  Sp  chosen  so  that  the  ground  state  energy  of  the  cylindrical  ion  coincides  with  the 
ground  state  energy  of  the  Coulombic  spherical  ion  with  the  same  Z.  The  product  UqA0 
is  chosen  so  that  the  corresponding  electric  held  coincides  with  the  barrier  suppression 
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Parameter 

Symbol 

Value 

Initial  State 

jz) 

|o,|> 

Steps 

Nt 

216 

Cells 

NX  X  Ny 

to 

X 

to 

Time  Step 

At  0.03 h/mc2 

Space  Step 

Ax,  Ay  or  Ax,  A z 

0.12  h/mc 

Residual  Charge 

Z 

54 

Soft  Core  Radius 

5p 

1.9  h/mc 

Laser  wavelength 

A 

1.52  nrn 

Vector  Potential 

^4o 

2.6  mc2/e 

TABLE  I:  Parameters  for  simulation  of  relativistic  ionization 

threshold  model  [5].  There  is  a  constraint  imposed  by  limited  computation  time,  which 
tends  to  increase  ujq.  For  the  example  given  here,  cuo  is  chosen  to  make  the  simulation 
time  a  few  hours  on  a  16-GPGPU  cluster,  while  keeping  the  adiabaticity  parameter  in  the 
tunneling  regime.  The  simulation  parameters  are  given  in  Table  I.  The  initial  bound  state 
is  computed  by  solving  (8)  and  resampling  the  result  on  the  Cartesian  grid  of  the  time 
dependent  problem. 

The  simulated  ionization  rate  is  defined  in  terms  of  the  charge  current  flowing  out  of  a 
volume  containing  the  bound  state.  Solutions  of  the  Dirac  equation  satisfy  the  continuity 
equation  <9/tj7t  =  0  where  j11  =  In  three  dimensional  form,  dtQ+  V  •  j  =  0,  where 

g  =  (32) 


The  expectation  value  of  the  charge  contained  in  a  ball  B  is 

(q)b  =  [  (33) 

Jb 

Take  the  radius  of  B  to  be  large  enough  so  that  it  contains  the  charge  associated  with  the 
bound  state  almost  entirely.  Then  the  ionization  probability  is  defined  as  the  expectation 
value  of  the  charge  outside  B,  divided  by  the  total  charge  of  the  particle: 


v  =  (g)s 
q 


(34) 
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FIG.  6:  Simulated  and  analytical  ionization  rates.  The  solid  curve  is  the  simulated  rate,  defined 
as  in  (35),  where  B  is  a  cylinder  of  radius  20A.  The  dashed  curve  is  based  on  the  dressed  Coulomb 
corrected  S'-matrix  analysis  of  Ref.  [6] .  In  plotting  the  dashed  curve  a  factor  was  inserted  to  unwind 
cycle  averaging. 


The  ionization  rate  is 


...  <ip  i  n  r ,, 

w  =  m  =  -qmjeiie 


(35) 


which  by  the  divergence  theorem,  is  the  same  as  the  current  flowing  out  of  the  volume, 
divided  by  the  charge. 

It  is  of  interest  to  compare  the  simulated  ionization  rate  with  various  analytical  formulas. 
In  this  report  we  make  comparison  with  the  dressed  Coulomb  corrected  S'-matrix  theory 
of  Ref.  [6].  Strictly,  the  S- matrix  theory  gives  the  ionization  rate  averaged  over  an  optical 
period.  However,  in  the  adiabatic  limit,  the  cycle  averaging  can  be  unwound,  and  the  rate 
can  be  considered  an  instantaneous  function  of  the  field  [7-10] .  The  result  of  this  procedure 
is  shown  in  Fig.  6.  The  simulated  rate  is  larger  than  that  of  the  S-matrix  prediction.  The 
peak  simulated  rate  is  delayed  with  respect  to  the  peak  of  the  electric  field  (93  =  0)  because 
of  the  propagation  delay  between  the  ionic  core  and  the  diagnostic  sphere  d B. 

A  useful  measure  of  the  accuracy  of  any  simulation  of  photoionization  is  the  precision 
with  which  charge  is  conserved.  Fig.  7  shows  the  evolution  of  the  charge  error,  normalized 
to  the  electronic  charge,  during  the  course  of  the  simulation.  The  accuracy  is  about  6  digits 
until  significant  ionization  sets  in.  The  accuracy  drops  to  about  4  digits  by  the  end  of  the 
simulation. 
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FIG.  7:  Accuracy  of  charge  conservation  during  the  course  of  the  simulation  of  relativistic  tunneling 
ionization.  The  error  in  the  total  charge  is  denoted  5Q. 

Finally,  a  visualization  of  the  ionizing  wavefunction  is  shown  in  Fig.  8.  The  electric 
held  points  in  the  — ^/-direction  and  the  radiation  wavevector  points  in  the  ^-direction.  As 
expected,  the  ionized  wavepacket  is  extracted  opposite  the  direction  of  the  electric  field,  and 
is  also  pushed  in  the  direction  of  the  radiation  momentum. 


IV.  COMPUTATIONAL  PERFORMANCE 

The  programming  model  used  for  the  Dirac  equation  follows,  almost  without  modification, 
the  methods  given  in  Ref.  [1]  for  the  Klcin-Gordon  equation.  The  primary  difference  is  that 
in  the  case  of  the  Dirac  equation,  additional  synchronization  is  necessary  to  guarantee  that 
the  electron  state,  0,  is  updated  before  advancing  the  positron  state,  x ■  The  simplest  way 
to  do  this  is  to  use  two  separate  OpenCL  kernels  to  advance  (p  and  %.  Apart  from  this 
distinction,  all  the  methods  for  distributing  the  workload  on  a  cluster  of  GPGPUs  remain 
the  same. 

The  floating  point  precision  used  in  the  code  can  be  selected  at  compile  time.  Single 
precision  is  suitable  for  cases  where  one  is  interested  in  a  final  state  which  is  significantly 
different  from  the  initial  state.  If  the  initial  state  is  weakly  perturbed,  accuracy  typically 
demands  that  double  precision  to  be  used. 
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FIG.  8:  Charge  density  associated  with  ionizing  relativistic  wavefunction  evaluated  at  (a)  mc2t/h  = 
490,  (b)  mc2t/h  =  980,  and  (c)  mc2t/h  =  1470. 

The  performance  on  various  GPGPU  devices  is  given  in  Fig.  9,  for  single  and  double 
precision,  using  a  2000  x  2000  numerical  grid.  The  performance  penalty  for  using  double 
precision  is  about  two-fold  on  the  gaming  device  (NVIDIA  Geforce  GTX  980),  and  somewhat 
less  on  the  other  devices.  The  water  cooled  device  advances  about  350  million  cells  per 
second.  Performance  could  likely  be  improved  by  optimizing  memory  access  patterns. 

The  performance  on  a  cluster  of  GPGPU  devices,  namely  the  SGI  ICE  X  system  “Topaz” 
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FIG.  9:  Performance  in  terms  of  millions  of  cells  per  second  for  various  GPGPUs.  Blue  bars  are 
for  single  precision  and  green  bars  are  for  double  precision.  The  Firepro  W9100  is  water  cooled. 
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FIG.  10:  Strong  scaling  study  on  the  DoD  HPCMP  cluster  “Topaz”  in  terms  of  millions  of  cells 
per  second  in  the  aggregate.  Double  precision  is  used  in  all  cases.  One  GPGPU  is  assigned  to  each 
MPI  node. 

at  the  U.S.  Army  Engineer  Research  and  Development  Center  (ERDC),  is  given  in  Fig.  10. 
The  Topaz  system  has  32  GPU  nodes,  each  with  a  single  NVIDIA  Tesla  K40P  device. 
Fig.  10  gives  scaling  results  for  a  2000  x  2000  numerical  grid.  In  going  from  16  to  32  nodes, 
performance  is  still  gained,  so  the  strong  scaling  limit  is  not  reached.  In  order  to  carry  out 
a  meaningful  weak  scaling  study,  a  larger  GPGPU  cluster  would  be  needed. 
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V.  CONCLUSIONS 

Numerical  solution  of  the  Dirac  equation  is  applied  to  the  problem  of  relativistic  tunneling 
ionization.  Analytical  and  numerical  solutions  of  the  time  independent  equation  are  used 
to  obtain  the  initial  wavefunction  as  a  bound  state  of  a  soft-core  potential.  An  efficient 
numerical  integration  scheme  for  the  time  dependent  equation  is  developed  to  advance  the 
initial  wavefunction  in  the  presence  of  a  relativistically  intense  electromagnetic  wave.  The 
numerical  scheme  is  implemented  to  take  advantage  of  clusters  of  GPGPU  devices.  It  is 
found  that  the  numerical  ionization  rate  exceeds  the  dressed  Coulomb  corrected  SFA  rate. 
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